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high energy neutrinos. The method is based on the maximisation of the detection power maintaining 
the control of the confidence level of an hypothetical discovery. This is achieved by using the so-called 
"False Discovery Rate" (FDR) controlling procedure. It has the advantage to be independent of the 
signal modelling and to naturally take into account the trial factor. Moreover it is well suited to the 
detection of multiple sources. 
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1. Introduction 

With the construction of the IceCube detector [1 
jin the ice of the South Pole glacier following 10 
years of data taking by its predecessor AMANDA^ 
and the completion of Antares [3] in the depth of 
the Mediterranean sea opening the way to its kilo- 
metric successor Km3NeT0, high energy neu- 
trino astronomy is entering a new era. 
As neutrinos are neutral and only sensitive to 
weak interaction they can directly be related to 
their production source. Their detection would 
be a proof of the existence of hadronic processes 
associated to the most violent objects of the Uni- 
verse and could provide the so far missing "smok- 
ing gun" for the acceleration of high energy cosmic 
ray protons. The neutrino telescopes are based on 
the detection of Cherenkov light emitted by the 
secondary lepton produced by the charged cur- 
rent interaction of neutrinos in a large volume 
(~ 1 km 3 ) of transparent medium (ice or sea 
water) using an array of photomultipliers. The 
present analysis focuses on the detection of muons 
produced by charged current interaction of high 



energy (> 100 GeV) neutrinos. The muon track 
direction is reconstructed by using the arrival time 
of Cherenkov photons and the muon energy is es- 
timated using the number of triggered photomul- 
tipliers. The neutrino direction accuracy depends 
on the kinematics of the reaction and the intrinsic 
detector resolution (of the order of 1° in the ice 
for AMANDA and 0.1° in the water for Antares). 
These detectors have to cope with a physical back- 
ground produced by the interactions of high en- 
ergy cosmic rays with the atmosphere. The prompt 
muons are suppressed by using the Earth as a 
shield while the atmospheric neutrinos constitute 
an irreducible background. The aim of the present 
analysis is to identify spatially localised excesses 
of events on top of this background. 
In this paper, an original approach is presented 
to increase the sensitivity to high energy neutrino 
point source signal, while controlling the confi- 
dence level of an hypothetical discovery in a model 
independent way. We will focus on time inte- 
grated northern sky surveys, looking for sources 
independently of known electromagnetic counter- 
parts. The current statistical approaches and their 
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limitations as well as the basics of the False Dis- 
covery Rate controlling procedure are presented in 
section 2. The application of FDR to the search 
of point sources with the AMANDA neutrino tele- 
scope is presented in section 3 and the perfor- 
mances of the method in section 4. 

2. Statistical tools for point source search 

2.1. The multiple hypothesis tests 

Looking for a point source in a sky of events 
can be seen as multiple hypothesis tests. Indeed, 
events which are described by some random vari- 
able X can either come from background and fol- 
low a statistical law Hq called the null hypoth- 
esis, or be signal events following an alternative 
Hi law. A testing procedure consists in defining 
a so-called critical region which is a sub-region of 
the possible values of X where the null hypothe- 
sis will be rejected and hence the related events 
considered as sources. Outside this region, events 
will be considered as background. The possible 
different outcomes are summarized in Table [TJ 



it will lead to over conservative tests. Another in- 
dicator, that we will use in the following, is the 
False Discovery Rate (FDR) defined as the expec- 
tation value of the rate of false rejections among 
all rejections: 



FDR = E( 
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A simple interpretation is that the expectation 
value of the confidence level of a discovery will 
be 1 — FDR, since FDR is an estimator of the 
size of the test. Another important quantity is 
the number of Type II errors S, the non-rejection 
of the signal. The probability of not making such 
an error at each test defines the power of the pro- 
cedure. One is generally interested in finding the 
testing procedure which maximizes the power for 
a given size. Technically, such a procedure will 
consist in finding a variable to test and defining 
the critical region where the null hypothesis will 
be rejected. In most cases the critical region will 
be characterized by setting an adequate threshold 
on this variable. 
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Table 1: Summary of the different outcome counts of mul- 
tiple hypothesis tests. 

Several "quality" indicators for a given test- 
ing procedure can be constructed from these num- 
bers, first by the number of Type I errors U which 
are the (false) rejections of true nulls. The expec- 
tation of U is called the size of the test. The 
most widely used related quality indicator of a 
testing procedure is the Family Wise Error Rate 
(FWER) which is defined as the probability of 
making at least one Type I error: 



FWER = Pr(U > 0) 



This indicator might appear inadequate in case a 
source is made up of several rejected events where 



2.2. Classical approaches and their limitations 

When both Hq and H\ are known, the most 
powerful procedure is the Neyman-Pearson test 
which relies on likelihood ratio Q. One defines: 
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where £ is the likelihood function for the N out- 
comes Xi of the tested random variable X consid- 
ering hypotheses Hj with j G {0, 1}. The critical 
region will be defined by A < c a where c a is such 
that if one desires a confidence level of at least 
1 — a: 



Pr(X(x) < c a \H ) 



a 



(4) 



This is what is commonly used in neutrino astron- 
omy [3] [19]. But there are important drawbacks. 
The first problem of this approach is its depen- 
dence on signal models since one has to know 
the Hi probability density function (pdf) mak- 
ing the method strongly model dependent. The 
second is that the determination of c a is often not 



explicitly computable and needs to be evaluated 
using heavy Monte Carlo simulations. Indeed in 
practice when one wants to insure a high con- 
fidence level the critical region will be character- 
ized by very low statistics. One usually overcomes 
this problem by extrapolating the distribution of 
higher p-values. Moreover, instead of looking for 
a realization below an a priori fixed threshold, 
one usually looks for the highest excess in data 
and then checks the probability to obtain a lower 
p- value, adapting this way the number of Monte- 
Carlo equivalent experiments to the data. The 
problem here is that it is not suited to the detec- 
tion of multiple sources. 

Alternatively a model independent "blind" anal- 
ysis jl] could be performed using only the back- 
ground hypothesis which can be inferred from the 
data. This method avoids the use of heavy Monte 
Carlo simulations for the determination of the 
critical region. Furthermore, fixing a priori the 
number of false nulls does not limit the multiple 
source detection capability. 

In the following we will use a set of N realizations 
{%i}ie{i,..,N} of a random variable X, and its pdf 
Po(x) under the null hypothesis H . The first step 
is to compute p-values p. t for each realization Xi, 
defined as the probability to observe at least Xi 
under the null hypothesis: 



Pi 



Pn(x) dx 



(5) 



A naive way to perform the test would be to 
reject all realizations with p-value lower than a 
fixed threshold pt = a. But after iV tests the fi- 
nal confidence level of any discovery will be lower 
than (1 — p t ) x N in the most general case, and 
1 — (1 — p t ) N in the particular case of indepen- 
dent tests. This is the trial factor effect which 
we would like to take into account without large 
numbers of heavy Monte-Carlo simulations. 
A very conservative way to overcome this and 
hence to control the FWER is the so-called Bon- 
feronni approach 10j. It consists in dividing the 



individual threshold by the number of trials (the 
equivalent for the independent case is the Sidak 



(1 — (1 — Pt) Jf )- But then the detection power 
decreases very strongly since the individual mini- 
mum threshold p-value will be divided by N, lead- 
ing to an increase of the number of type II errors. 
We will describe in the following a procedure which 
controls the FDR and hence the confidence level, 
and which by relaxing the FWER control makes 
it possible to reach a high detection power. 

2.3. False Discovery Rate controlling procedure 
This FDR procedure has been first developed 



11] and applied to several kinds of astrophys- 
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ical and cosmological searches 
based on an adaptive threshold on p-values which 
guarantees both high detection power and con- 
trol of the FDR and hence the confidence level of 
an hypothetical detection. The procedure itself is 
rather simple. Given a set of N p-values and a 
FDR input value a: 

• sort the p-values in increasing order, to get 
the ordered set {pi}i^i,..,N- 



find i c so that 



1 i 

i c = max(i e {1, .., N} | p < < -.a.— )(6) 

X N 



where x is a coefficient that has to be intro- 
duced to account for dependency between 
the tests. It is generally 1 and the the- 
oretical value in the worst case is J2iLi \- 
The latter value makes the procedure more 
conservative but is rarely necessary. In the 
particular case of p-values distributed uni- 



formly the coefficient is unity 15 



• all tests with index i < i c will reject the Null 
Hypothesis and we considered as sources. 

The procedure is proved to be "minimax" (op- 
timal choice minimizing risk in decision theory) 
for Gaussian, Poissonian and exponential distri- 



butions [16] . Moreover, unlike the classical ap- 
proach it is well suited for the detection of mul- 
tiple sources since it relies on the rate of false 
discoveries. 



10] procedure which replaces the threshold by 
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3. Application to AMANDA reconstructed 
high energy neutrino skies 

The FDR method described in what follows is 
based on the use of a random variable built as a 
convolution of 2 variables: the angular position 
and the energy of the neutrino candidates. It has 
been tested on event distributions corresponding 
to the high energy neutrino northern sky as seen 



signal ratio r^, s , corresponding [20[ to 1.585 x p 
. If fimax is the maximum value of p>bg{$) over the 
sky, in order to ensure the same background event 
count n m ax all over the sky R must be: 



R(S) = Arccos(l- 



'l-cosR opt (S)))(t 



by AMANDA [18f] . These data have been used to 
determine the atmospheric neutrino background 
sample needed to construct the Null Hypothesis 
pdf. On the other hand a cosmic neutrino sample 
was used to derive limits and discovery potential 
corresponding to one year of data taking as well 
as to compare with the classical method based on 
Likelihood ratio. 

3.1. Null Hypothesis 

3.1.1. Information from angular position 

Due to the detector cylindrical symmetry along 
the north-south axis, the distribution of the an- 
gular position is uniform in right ascension and 
therefore only depends on the declination 5. Given 
the low statistics accumulated annually by the 
AMANDA detector (of the order of 1000 events), 
the number of neutrinos in a sky area of the order 
of the detector resolution p ( ~ 3 deg ) follows a 
Poisson statistic. For each reconstructed neutrino 
in this area corresponds a circular search region 
for which the space angle radius R will vary as 
a function of 5 in order to ensure the same Null 
Hypothesis all over the sky. 

The detector resolution is defined as the median 
of the point spread function (PSF) which is the 
distribution of the differences between the true 
neutrino direction and the reconstructed one. Ne- 
glecting small dependency on the neutrino energy, 
the angular resolution is assumed to be the same 
for signal and background events. The mean num- 
ber of expected background events in a region of 
angular radius R opt (S) is defined as: 



This way R will be at least as large as R op t- This 
is the best compromise since for a circular area 
the ratio r^i s increases slower for values higher 
than R opt than for lower ones as can be seen in Fig 
[1] which represents the evolution of r^, s with the 
radius for an hypothetical PSF of 1° resolution. 




Evolution of f^, g versus the radius of the 



Figure 1 

bin for a resolution of 1 



3.1.2. Information on energy 

The energy is another observable discriminat- 
ing signal and background, since the cosmic spec- 
trum is expected to be generically harder than the 
atmospheric one. Cosmic neutrino energies are 
supposed to follow a power law of index —2 ver- 
sus —3.7 for atmospheric secondaries. The most 
robust and simple variable available to evaluate 
the neutrino energy is the number of hits N c h, de- 
fined as the number of optical modules triggered 
by the photons emitted by the muon when pass- 
ing through the detector |6| . Like the background 
angular position distribution, the pdf Pq(N c}i ) cor- 
ing to this variable can be extracted from 
the data (Figure [3]). It is combined with the angu- 
where R opt (5) is the angular radius of the region lar position distribution to lead to the total energy 
minimizing the square root of background over accumulated in a search region of the sky. 
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1 - cos R opt (5) 
2 cos 5 sin R opt (5) 



S+R opt (S) 



8-R opt (S) 



In the case of a single neutrino the only informa- 
tion is its own energy which pdf is Pq. In the case 
where there is only one neighbour neutrino, the 
total energy distribution will be : 



(9) 



a=0 



By recursion, the total accumulated energy within 
a search region centered on a neutrino and con- 
taining n extra neutrinos is distributed as : 



P n (x) = J2 P o(a) P n -i{x -a) 



(10) 



a=0 



These probabilities have to be weighted by the 
Poissonian distribution V{n\ fJ, ma x) to observe n 
neutrinos in a region where fi ma x are expected. 
The total energy distribution corresponding to 
the Null Hypothesis is: 

oo 

PTOT{N ch ) = V ( U ' Vmax) P n {N ch ) (11) 
n=0 

3.2. Clustering 

The method can lead to multiple counting of 
signal events after the application of the FDR pro- 
cedure as the region attached to each neutrino can 
overlap. One way to overcome this is to perform 
a spatial clustering on rejected neutrinos. The 
most suitable scheme here is a hierarchical min- 
imal distance clustering 21] . The algorithm can 



be summarized as follows: 

Initially each neutrino of the rejected set N r is at- 
tached to a different cluster Q. 
The following procedure is then applied as long 
as the number of clusters is decreasing. 
If the distance between two clusters is lower than 
the sum of the resolutions R(a) and R(b) on the 
direction of the corresponding neutrinos, the two 
clusters they belong to are merged 

if dij < R(a) + R(b)then d = Cj = Q U Cj 
where 

dij = {mm(distance(v a , Vb)) I v a £ C« , G Cj} 

At the end of the clustering, each cluster will 
be considered as a signal source. The false discov- 
ery rate estimator does not account for individual 
neutrinos anymore but for clusters. 



3.3. Final cut 

The procedure at this stage does not yet al- 
low perfect control of the false discovery rate at 
high source luminosity. This can be understood 
from the principle of the FDR rejection, which is 
designed to detect tiny excesses with respect to 
the Null Hypothesis. Hence a source with high 
potential of detection obviously deviates from the 
main purpose of the method. This generates arti- 
facts appearing as small clusters mainly made up 
of a single neutrino. For consistency, we want the 
method to control the False Discovery Rate even 
in these extreme cases. The simplest and most 
efficient way is to cut off single neutrino clusters 
when at least another bigger cluster exists within 
the same events sample. This does not affect the 
efficiency of the method at low luminosity. 

4. Results 

The method has been tested on simulation 
samples where neutrino point sources have been 
added to an atmospheric background of 1000 events 
corresponding to the statistics collected by AMANDA 
in one year 13]. The statistical estimators of false 
discovery rate and detection power are defined as 
the mean value of each quantity over a large num- 
ber of artificial sky realizations. Given the bino- 
mial nature of these estimators (a discovery is ei- 
ther true or false) the statistical error e will be 
H: 



x(x — 1) 
N — 1 



(12) 



where x is the estimator and N the number of sky 
realizations used to determine the mean. With 
10 3 sky realizations one reaches a statistical er- 
ror of order 0.22% which is sufficient to test a 
99.5% confidence level control. The following sec- 
tion first describes the variables and pdfs used 
to define the Null Hypothesis and to build back- 
ground and source Monte Carlo generator. The 
false discovery rate control by the procedure is 
then checked and discovery potential of the method 
is quantified. 
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4-1. Generation of the background test sample 

As described in section 13. 1[ the two observ- 
ables used by the method are the angular position 
of events and their energy. Hence the distribu- 
tions of the declination and of the energy of the 
real events were fitted with analytical curves. The 
expected background as a function of the sine of 
the declination is shown in Figure [2] It has been 
fitted with a polynomial function. 

The number of hits (N c h) triggered by the pass- 
ing muon is used as energy estimator. Its distri- 
bution is shown in Figure [3] together with the fit 
used for the Monte-Carlo generation. 




sin (declination) 



Figure 2: PDF of the sine of the declination of events 
observed with the AMANDA detector. 

Since the detector is symmetric around the 
north-south axis, the right ascension distribution 
is uniform. These three pdfs are used to generate 
the Monte Carlo realizations of background skies. 

4-2. Generation of the test source sample 

Neutrino sources were generated on the basis 
of the detector response corresponding to an en- 
ergy distribution following a power law of index 
—2. These events were submitted to the full re- 
construction chain of AMANDA. The energy dis- 
tribution of these events is represented in Fig. |3j 
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Number of Hits 

Figure 3: Normalized distributions of the number of 
hits for the AMANDA detector. Points are data 
(and hence considered as purely atmospheric) and the 
line is the fit of the corresponding probability density 
function. Empty squares correspond to a simulation 
of the detector response to a cosmic flux. 

For all tests, one source is simulated with lumi- 
nosity varying between (pure background) and 
20 neutrinos, at various declinations to be com- 
pared with an average expected background count 
of two. In order to take into account the detec- 
tor resolution, cosmic neutrinos are spread around 
the source location following the PSF at this dec- 
lination. 

An example of point spread function can been 
seen in Figure H] and the evolution of its charac- 
teristic width in Figure E] ■ 

4-3. Construction of the Null Hypothesis pdf 

As mentioned in section [3XTJ a circular search 
region of angular radius R opt is centered on each 
reconstructed neutrino. The expected number of 
events (Figj5j) is computed thanks to the known 
background density. The search region radius R(5) 
is then defined in such a way that it contains al- 
ways the maximal expected background of 1.72 
neutrino independent of its angular position. The 
variation of this radius as a function of the decli- 
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Figure 4: Point spread function at 22.5° of declina- 
tion. The curve shows the fit function used for MC 
generation. 
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Figure 5: Expected background in optimal bin size 
Ropt- Maximal value is used as expected background 
to build bin size in the whole sky. 

nation can be seen in Figure [6j The pdf of the 
Null Hypothesis is then computed from EqJTT] 
The corresponding cumulative density function is 



Figure 6: Detector resolution R , optimal signal to 
noise ratio bin size R op t and the optimal constant 
background bin size S which is eventually used as 
bin definition. 

shown in Figure [71 




Figure 7: Cumulative Distribution Function on Total 
number of hits expected for a pure background sky 
with 1000 atmospheric neutrinos 
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44- FD R control 

The false discovery rate control was checked at 
ten different source declinations varying between 
0° and 85° and with source luminosities ranging 
from to 20 neutrinos. The false discovery rate 
is well controlled in all configurations at the re- 
quested 99.5% confidence level. An example can 
be seen in Figure [8] for a source located at 22.5° 
of declination. 
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Figure 9: Detection probability of a source located at 
22.5° above a background of 1000 atmospheric neutri- 
nos versus source luminosity. Grey curve corresponds 
to discrete number of detected neutrinos while the 
black one corresponds to a Poisson averaged mean 
luminosity. 

n neutrinos weighted by the Poisson probability 
V(n\h> STC ). 



Figure 8: False Discovery Rate obtained with 1000 
atmospheric neutrinos versus luminosity of a source 
located at 22.5° declination and requested FDR of 
0.5%. 



4-5. Discovery potential 

The probability to identify the source is esti- 
mated on the same test sample. Results are shown 
in Figure for a source located at 22.5° of dec- 
lination. The grey curve corresponds to the de- 
tection probability VV ev (u src ) of a given number 
of detected neutrinos from the source v STC . The 
number of detected events emitted by a source 
of a given luminosity follows a Poisson distribu- 
tion with a mean depending on the source flux. 
Consequently the detection power (black curve in 
FigJHD function of the luminosity u STC is com- 
puted by summing the probabilities of detecting 
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(13) 



To get physically meaningful information these 
numbers of source neutrinos have to be translated 
into fluxes by performing a convolution with the 
effective area of the detector. The fluxes needed 
to have 50% (FigM) and 90% (FigEJ chance to 
detect a source with a confidence level of 99.5% for 
one year of data-taking were computed. For in- 
stance, as can be seen in FigJHl the minimal mean 
source strength is 7.1 cosmic neutrinos over 1.72 
background neutrinos to have 50% probability to 
identify the source at a declination of 22.5°. This 
number increases to 12.6 cosmic neutrinos for a 
detection probability of 90%. 

The comparison of these results with other 
methods is not straightforward since, the FDR is 
able to handle multiple sources while other meth- 
ods, eg. likelihood ratio maximisation are opti- 
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mized for a flux model and only deal with the 
highest excess on a sky. However, it is reasonable 
to compare log likelihood ratio (LLH) method cor- 
responding to a 5a excess above background and 
FDR procedure at 99.5%CL. Indeed the trial fac- 
tor effect for a statistic of 1000 neutrinos in a sky 
will give a probability of approximately 0.5% that 
an excess corresponding to at least 5a would ap- 
pear on a pure background sky [23j | . 
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Figure 10: Flux discovery potential for a detection 
efficiency of 50% after 1 year of data taking 

For comparison, fluxes corresponding to de- 
tection efficiencies of 50% and 90% are presented 
in Fig {10] and [TTl together with the corresponding 
curves for LLH for 3a and 5a excesses. It has to 
be underlined that due to the trial factor effect an 
excess of 3a or higher would appear with a prob- 
ability of 64% in a pure background sky of 1000 
events. It can be seen that the FDR method is 
clearly more sensitive than the LLH for a 5a ef- 
fect at low declinations. This behaviour is linked 
to the increase of background in this region of the 
sky which implies for the LLH method an increase 
of the source luminosity to maintain the same de- 
tection efficiency 
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Figure 11: Flux discovery potential for a detection 
efficiency of 90% after 1 year of data taking 



4-6. Limits and sensitivity 

We use the unified ordering algorithm pro- 
posed by Feldman and Cousins [24j to construct 
confidence belts and quantify the sensitivity of 
the FDR procedure. Sensitivity is defined here 
as the average upper limit obtained from confi- 
dence belts for all possible experimental outcomes 
weighted by their probability if no signal is present. 
An example of a confidence belt is shown for a 
possible source located at 22.5° declination in Fig- 
ure [T2]. The sensitivity shown in this figure cor- 
responds to a differential cosmic neutrino flux of 
0.77 x 10 -10 TeV cm~ 2 s _1 . The variation of the 
sensitivity as a function of the declination is pre- 
sented in Figure [13j 

5. Summary and conclusions 

We have presented an original statistical ap- 
proach to the search of point-like sources of cos- 
mic neutrinos with Cherenkov telescopes, which 
allows to control the confidence level of an hypo- 
thetical discovery with multiple hypothesis test. 
Its main advantages are that it is independent of 
any model of the signal production, that it takes 
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Figure 12: Feldman- Cousins confidence belt for decli- 
nation 22.5° after 1 year of data taking 
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Figure 13: Sensitivity and discovery potential after 1 
year of data taking 

into account the trial factor effect and that, un- 
like existing methods, it naturally includes the 
search for multiple sources. Its performances were 
compared with a likelihood ratio minimisation ap- 
proach using the data collected by the AMANDA 



neutrino telescope. In the case of the detection 
of a single source in the context of a particular 
flux model, the FDR method results in a sensitiv- 
ity equivalent to the Bayesian approach, showing 
even slightly better performances at low declina- 
tion or in case of a very faint source. 
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